Ising nematic phase in ultra-thin magnetic films: a Monte Carlo study 
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We study the critical properties of a two-dimensional Ising model with competing ferromagnetic 
exchange and dipolar interactions, which models an ultra-thin magnetic film with high out-of-plane 
anisotropy in the monolayer limit. We present numerical evidence showing that two different sce- 
narios appear in the model for different values of the exchange to dipolar intensities ratio, namely, 
a single first order stripe - tetragonal phase transition or two phase transitions at different temper- 
atures with an intermediate Ising nematic phase between the stripe and the tetragonal ones. Our 
results are very similar to those predicted by Abanov et al [Phys. Rev. B 51, 1023 (1995)], but 
suggest a much more complex critical behavior than the predicted by those authors for both the 
stripe-nematic and the nematic-tetragonal phase transitions. 

We also show that the presence of diverging free energy barriers at the stripe-nematic transition 
makes possible to obtain by slow cooling a metastable supercooled nematic state down to tempera- 
tures well below the transition one. 

PACS numbers: 75.70.Kw,75.40.Mg,75.40.Cx 
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I. INTRODUCTION 

The study of thin magnetic films have deserved an in- 
creasing interest during the last decade, both on their 
experimental and theoretical aspects. Besides the great 
amount of technological applications related to their 
magnetic behavior, as for instance, data storage, these 
studies have also faced statistical physicists with the chal- 
lenge of trying to answer many foundational question re- 
garding the role of microscopic competing interactions 
on the macroscopic static and dynamical behavior of two 
dimensional systems. Moreover, the constant develop- 
ment of novel methods for constructing ultra thin films, 
together with a significant improvement in the quality 
of techniques for measuring nanoscopic magnetic struc- 
tures, e have open many fascinating questions, many of 
them still open. 

Many ultrathin magnetic films, Hke e.g. Co/Cu, 
Co/Au, Fe/Cu, undergo a reorientation transition at a 
temperature Tr in which the spins align preferentially 
in a direction perpendicular to the film iSii. This re- 
orientation transition is due to the competition between 
the in-plane part of the dipolar interaction and the sur- 
face anisotropyi. Furthermore, in the range of tem- 
peratures where the magnetization points out of the 
plane, the competition between exchange and dipolar 
interaction causes the global magnetization to be effec- 
tively zero but instead striped magnetic domain patterns 
emerge^. In the limit when the stripe width is much 
larger than the domain walls, the walls can be approx- 
imated by Ising walls and the system can be consid- 
ered as an Ising system of interacting domain walls^. 
In spite of intense theoreticali?i'^i''?i^iTifji'?|ifti"|i^|ii?|i^ and 
experimentali*2iiSiiLii work on the behavior of ultrathin 



magnetic films, the precise nature of the phases and the 
relaxational dynamics aspects of these systems is still 
poorly understood. In Monte Carlo simulations of an 
Ising model on a square lattice. Booth et alffi> found ev- 
idence of a stripe phase a low temperatures, with orien- 
tational and positional order reminiscent of the smectic 
order in liquid crystals. These authors found a transi- 
tion from a stripe phase to a high temperature phase 
with broken orientational order, called tetragonal liq- 
uid phase. In this phase domains of stripes with mu- 
tually perpendicular orientations emerge and form kind 
of labyrinthine patterns. At still higher temperatures 
these domains collapse and the system crosses over con- 
tinuously to a completely disordered behavior. From nu- 
merical data for the specific heat these authors found 
evidence of a continuous stripe/tetragonal-Hquid transi- 
tion. However, it has been argued that the transition 
observed by Booth et al should correspond to a nematic- 
tetragonal phase transitioni. Moreover, recent numeri- 
cal results using time series methods on the same model 
with the same parameters values suggest that such tran- 
sition is not continuous but it is a first order onei^. In 
an important contribution Abanov et ali^ analyzed the 
different phases which could emerge from a continuous 
approximation to the domain wall crystal. They pre- 
dicted two possible scenarios in zero magnetic field. In 
the first scenario a smectic-like low temperature phase 
with spatial correlations decaying algebraically with dis- 
tance appears at low temperatures. In this phase there 
is quasi long range order (QLRO) and the orientational 
order parameter is non zero. The proliferation of bound 
dislocations pairs at higher temperature causes QLRO to 
be destroyed through a Kosterlitz-Thouless phase tran- 
sition to a nematic-Hke phase. In this phase only orien- 
tational order is observed. At still higher temperatures 
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even orientational order is destroyed through the appear- 
ance of domains of perpendicular stripes. This is the 
tetragonal-Hquid phase. Since Z2 orientational symme- 
try is restored at this transition the authors speculated 
that this transition should be in the Ising universality 
class. In the second scenario the system is not able to 
sustain a purely nematic phase and goes from the smec- 
tic directly to the tetragonal liquid through a first order 
transition. Recently Cannas et al.^^ found numerical ev- 
idence of this second scenario, also supported by a con- 
tinuous approximation with a Landau-Ginzburg Hamil- 
tonian, from which a fiuctuation-induced first order tran- 
sition is predicted. Up to now no evidence was found of 
the first scenario predicted by Abanov et al. In this paper 
we show that for large enough system sizes, Monte Carlo 
simulations of the same model studied in5J^> support the 
appearance of a sequence of two phase transitions of the 
kind predicted by Abanov et al. when the width of the 
equilibrium stripes is large enough, while for thin stripes 
the second scenario with a unique transition is found. 

We consider a system of magnetic dipoles on a square 
lattice in which the magnetic moments are oriented per- 
pendicular to the plane of the lattice, with both nearest- 
neighbor ferromagnetic exchange interactions and long- 
range dipole-dipole interactions between moments. The 
thermodynamics of this system is ruled by the dimen- 
sionless Hamiltonian^^: 

where 5 stands for the ratio between the exchange Jo > 
and the dipolar > interactions parameters, i.e., 
(5 — Ji^l Jd- The first sum runs over all pairs of near- 
est neighbor spins and the second one over all distinct 
pairs of spins of the lattice; ry is the distance, measured 
in crystal units, between sites i and j. The energy is 
measured in units of Jd- The overall (known) features 
of the equihbrium phase diagram of this system can be 
found in Refsi^i^iii*i^ii^iiS, while several dynamical prop- 
erties at low temperatures can be found in Refsi 7ii?i^fli^^ . 
The threshold for the appearance of the stripe phase in 
this model is 5c = 0.42 5t8.19 , As 5 increases the sys- 
tem presents a sequence of striped ground states, char- 
acterized by a constant width h, whose value increases 
exponentially with S^^. 

We show that energy and orientational order param- 
eter histograms present a sequence of three peaks for 
5 = 2, corresponding to three different thermodynamic 
phases, while for i5 = 1 only two peaks are observed, 
consistent with the presence of only one phase transi- 
tion. In the first case we show that the intermediate 
phase presents long range orientational order but no long 
range positional order, consistently with a nematic phase. 
Our results show that the low energy phase is a striped 
one. We did not find evidence of smectic order, although 
the possible existence of algebraic decaying correlations 
(near and below the transition temperature), strongly 
hidden by finite size effects, cannot be excluded. Finite 



size scaling analysis of specific heat and the fourth or- 
der cumulant of the energy give evidence of a first order 
transition between the nematic and tetragonal phases. 
However, the analysis of the orientational order param- 
eter histograms suggests the existence of more than one 
nematic phases for much larger system sizes than the 
ones considered here. On the other hand, the stripe- 
nematic phase transition shows unusual features, some of 
them characteristic of a first order transition, but some 
other properties strongly resembhng those observed in a 
Kosterlitz-Thouless (KT) phase transition. 

For i5 = 1 a unique weakly first order phase transition 
is supported by direct thermodynamic analysis through 
a computation of the free energy of the different phases. 
In the last section we show that the previous thermo- 
dynamic behavior is also supported by out of equilib- 
rium measurements during quasi-static cooling/heating 
cycles, where strong metastability is observed at the 
stripe-nematic phase transition. This behavior is con- 
sistent with the observation of asymptotically divergent 
free energy barriers at this transition. 

II. THE ISING NEMATIC PHASE 

We first analyzed the equilibrium histograms P{E/N) 
for the energy per spin at different temperatures T and 
different system sizes L. For every system size and ev- 
ery temperature the corresponding energy histogram was 
calculated by recording the energy values during a single 
MC run. Before starting to record the energy we left the 
system run a transient period of Mi MC steps (MCS) in 
order to equilibrate. After that period we recorded the 
energy values over M2 MCS.AMCSis defined as a com- 
plete cycle of spin update trials, according to a heat 
bath dynamics algorithm. For every pair of values of T 
and L we checked different values of Mi and M2 in order 
to ensure a stationary distribution P{E/N). Typical val- 
ues of Ml were between 10^ and 2 x lO'^ MCS, while typ- 
ical values of M2 were between 2 x 10^ and 2 x 10* MCS. 
A similar calculation was carried out to obtain the equi- 
librium histograms of the orientational order parametes^ 



ry = — (2) 

where n„ (nh) is the number of vertical (horizontal) 
bonds between nearest neighbor anti-ahgned spins. This 
parameter takes the value -1-1 (—1) in a completely or- 
dered horizontal (vertical) striped state, while it equals 
zero in any phase with 90° rotational symmetry, thus 
describing the 90° rotational symmetry breaking. 

We concentrated first in analyzing the 6 = 2 case {h = 
2 striped ground state). In FigQlwe show the behavior 
of P{E/N) for L = 56, 6 = 2 and different temperatures, 
while in FiglJlwe show the corresponding histograms for 
the distribution of the absolute value of the orientational 
order parameter P(|7y|). 
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FIG. 1: Energy per spin probability distributions (normalized 
histograms) for 5 = 2 and L — 56. 



The presence of pairs of peaks in both the energy per 
spin and the orientational order parameter distributions, 
located at three distinct ranges of values, is a clear signa- 
ture of the existence of three different phases, with two 
phase transitions between them. The low energy peak 
is associated with a peak in the order parameter cen- 
tered around \rj\ « 1, thus corresponding to a phase with 
long range orientational order. Moreover, an analysis of 
a stripe order parameter introduced below (and also a di- 
rect inspection of the typical spin configurations at those 
energies) shows that this is an ordered striped phase with 
h = 2. The highest energy peak is associated with a peak 
in the order parameter centered at 77 = 0. A direct in- 
spection of the associated typical spin configurations at 
those energies shows indeed that they correspond to a 
tetragonal liquid phase. On the other hand, the inter- 
mediate energy phase is associated with a peak of the 
order parameter centered around |ry| « 0.81, thus corre- 
sponding to a state with broken 90° rotational symme- 
try. In FiglSK we illustrate a typical spin configuration 
in this intermediate phase for L = 64 and T = 0.78. We 
see that this phase is characterized by a high density of 
topological defects, mainly dislocations in the directions 



12 - 

10 - 


T = 0.775 

1 




T = 0.780 

3 


^ 8 - 








6 - 


i 






4 - 








2 - 






A 


- 






12 - 
10 - 


T = 0.805 




T = 0.808 


^ 8 - 








6 - 

Q- 








4 - 








2 - 








- 








12 - 
10 - 


T = 0.815 




T = 0.830 


„ 8 - 








ar 6- 








4 - 








2 - 








- 









0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



hi |T1| 

FIG. 2: Orientational order parameter probability distribu- 
tions (normalized histograms) for 6 = 2 and L = 56. 

of the underlying striped structure. The same type of 
defects have been observed in Fe on Cu ultrathin films 
near the phase transition (see Fig. 2 in Ref?^). Such de- 
fects reduce the average value of the orientational order 
parameter and their presence is in agreement with the 
qualitative description of the Ising nematic phase intro- 
duced by Abanov and coauthorA To confirm this as- 
sumption we calculated the spatial correlations along the 
coordinate directions given by 

y X 

^f^*^) = -^^'^'\'^^,v^'^,v^r) (4) 

y X 

An example of the behavior of Cx {r) and Cy (r) in the 
intermediate phase is shown in Fig^ for L = 92 and 
T = 0.78. This figure shows clearly that the correlations 
decay algebraically in one of the coordinate directions 
and exponentially in the other. We see that the intermedi- 
ate phase is characterized by long range orientational or- 
der but does no t present positional order. Although the 
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FIG. 3: Typical spin configurations in the nematic phases 
for 5 = 2 and L = 64. (a) T = 0.78, the orientational order 
parameter for this configuration is |?7| ~ 0.8; (b) T = 0.807 
and « 0.59 
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form of (Tf7. It can be shown that in a pure striped 
state of width h the only non-zero components corre- 
spond to (kxjky) — {0,±7Tp/h) (horizontal stripes) or 
(fca;, ky) — (zLirp/h, 0) (vertical stripes), where p takes all 
the odd integer values between one and h. For instance, 
in a pure h — 2 vertical striped state it can be shown 
that 

'S'(fc) = ^ 4„,0 (4,,7r/2 + 4^,-7r/2) (6) 

Now, since all the phases we are dealing with in this 
case present a discrete rotational symmetry, either of 
90° or 180° respect to the coordinate axes, the maxima 
of S{k) will be located at {kx,ky) = (±fco,0) and/or 
{kx,ky) = (0,±A:o), with /cq some value close to it/2. 
Then we can define the stripe order parameter as 




FIG. 4: Spatial correlation functions along the coordinate 
directions in a nematic phase for 5 = 2, T = 0.78 and L — 92. 
The continuous line corresponds to a fitting using a func- 
tion /(r) = Aexp (r/^) sin (fco r + (f>), with fitting parameters 
ko = 1.47, C = 17.4 and (j> = 1.63. Notice that ko ~ n/2, 
the wave vector of the striped structure with h — 2 (ground 
state for 5 — 2). The log-log plot of Cx in the inset shows 
that it decays at large distances with a power law r~" , with 
an exponent uo ~ 0.12. 



description level of the elastic approximation of Abanov 
and coauthors does not allow a simple derivation of the 
spin-spin spatial correlations, the above features are in 
qualitative agreement with their characterization of the 
nematic phase. 

To further characterize the different phases we defined 
a stripe order parameter through the static structure fac- 
tor 



where cr^ = T^^p^^ '^f' is the discrete Fourier trans- 



77, = 2(^(fco,0) + 5(0,A:o)) 



(7) 



(this definition can be easily generalized to consider un- 
derlying ground states of arbitrary width h) . This order 
parameter will take values close to one for states with 
long range h — 2 stripe order and zero (in the thermo- 
dynamic limit) for any state without long range posi- 
tional order. For instance, let us consider a system in 
which correlations decay exponentially in the x direction 
and remain almost constant in the y direction, namely, 
C(r, r ') = exp (— A|a; — x'\) cos (A;o x); this is a rough ap- 
proximation of the behavior we observe in FigQJ where 
the algebraic decay in one of the coordinate directions is 
extremely slow. Substituting into the definition of S{k) 
and replacing the sums by integrals, in the large L limit 
we get 



S{k) ~ J^yj^ j g ■'^^cos{kQ x) coskx X dx 
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2L \ik.,-koy + X^ (fc^ + fco)' + A2 



8) 



If the correlation length — A^^ is independent of the 
system size we have that the stripe order parameter goes 
to zero as tJs when L — > oo. In a tetragonal liquid 

state we can assume for the correlation the following form 



C{r,f") = e~^(\^~^'\+\y-y'\') cos(ko \x-x'\) cos{ko \y-y'\) 

(9) 

which possesses the tetragonal symmetry. We have then 



S{k) 



A 



A 



1 

4^2 V (fc.. - fco)2 + A2 ' (fc, + fco)2 + A2 
A A 

{ky - fco)2 + A2 + {ky + fco)2 + A2 



(10) 



which reproduces the observed crown shape of the struc- 
ture factor observed in a tetragonal liquid state, both in 
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FIG. 5: (a) Static structure factors S{k,0) and S{0,k) at 
T = 0.78 for different system sizes. The dashed lines corre- 
spond to Lorentzian fittings f{k) = A/((fc — fco)^ + A^). The 
upper inset shows the fitted parameter kg vs. . The lower 
inset shows that the maximum of S{k) scales as . (b) 
2 L{S{k, 0) + 5(0, k)) vs k/ko for the same numerical data in 
(a); the continuous line is a Lorentzian fitting. The data col- 
lapse of the different curves shows that the whole curve S{k) 
scales as L^^, as predicted by Eq.ljSj. 



numerical simulations^ and in Fe/Cu films imagesi^. In 
this case we have that the stripe order parameter goes to 
zero as rjs ~ when L oo. 

We performed numerical calculations of the structure 
factor in the three different phases for different sys- 
tem sizes up to L = 72. First we calculated it in 
the low energy phase by letting the system to equili- 
brate from a vertical h = 2 ground state configuration 
at T = 0.77, just below the estimated transition tem- 
perature between the low and the intermediate energy 
phases Ti « 0.772 (see next section). We found that 

S{k) = c Sky,o (4,,7r/2 + 4,,-7r/2), whcrc c quickly con- 
verges to a value c « 0.48 consistently with Eq.ljSl, thus 
confirming the long range stripe character of the low en- 
ergy phase. Next we let the system to equilibrate in the 
intermediate phase at T = 0.78. In FigEt we plot S{k, 0) 
and S{0,k) for different system sizes. We observe that. 



FIG. 6: Static structure factors S{k, 0) and S{0, fc) at T = 
0.95 for different system sizes. The continuous lines corre- 
spond to Lorentzian fittings f{k) = X/{{k — ko)^ + A^). The 
upper inset shows the fitted parameter ko vs. . The lower 
inset shows that the maximum of S{k) scales as . It can 
be also shown that the whole curve S{k) scales as L , as 
predicted by Eq.JT^J. 



while one of the two functions remains almost zero for all 
k, the other presents two symmetric peaks at fc = i/co 
(only one is shown in the figure) which are very well fit- 
ted by a Lorentz function, in agreement with Eq.®. We 
see that kg saturates into a value ko ~ 1.48 while the 
maxima (or equivalently the order parameter rjs) goes to 
zero as L~^. Moreover, the data collapse of the curves 
for different values of L when 2 L{S{k,0) + S{0,k)) is 
plotted vs k/ko (see FigEJa) shows that the whole curves 
scale as in agreement with Eq.®; this confirms that 
correlation length does not depend on L in the thermo- 
dynamic limit and that the intermediate phase does not 
present LRO, thus confirming its nematic character. Fi- 
nally, in FigElwe show S'(fc,0) and S{0,k) for different 
system sizes after equilibration in the tetragonal Hquid 
at r = 0.95. We see that both functions present two 
symmetric peaks at A: = ±fco (only those at A: > are 
shown in the figure) which are very well fitted by Lorentz 
functions and that the maxima scale as in excellent 
agreement with Ea. Hl()|l . Hence, we see that the stripe or- 
der parameter allows, not only to distinguish long range 
orientational order, but also discriminate between the 
tetragonal and the nematic phases through its finite size 
scaling. 

It is worth mentioning that the two-peak structure in 
P{E) and P{r]) associated with the nematic-tetragonal 
phase transition can only be appreciated in both type 
of histograms only for system sizes L > 40. For system 
sizes L < 40 the two peaks are so close to each other 
that the phase transition cannot longer be resolved. In 
that situation the tetragonal and the nematic structures 
are almost indistinguishable, the latter appearing as a 
slightly elongated tetragonal onei^. While this strong fi- 
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nite size effect disappears for system sizes L > 40, there 
still remain similar finite size effects associated with the 
nematic-tetragonal phase transition up to system sizes 
L — 64. This can be seen in FigQ where we show the 
histograms for the energy and the order parameter for 
L = 64 and different temperatures around that transi- 
tion. If we look at the energy per spin histogram (FigE^) 
we just see that the high energy peak becomes skewed to- 
wards the right and broadens around T = 0.81 (compare 
with FigQJ. However, when looking at the orientational 
order parameter histogram (FigdD) we see that the sin- 
gle peak observed for smaller system sizes (see FigEJ at 
that temperature range splits for L = 64 into two dis- 
tinct peaks, one centered at « 0.8 and the other at 
I77I « 0.6. In FigEb we can see an example of a typical 
spin configuration in the second case. We see that the 
system still exhibits nematic order, the main difference 
with the configuration of FigEt being a higher density 
of domain walls perpendicular to the underlying striped 
structure (disclinations), which reduce the value of \ri\. 
These results suggest the existence of different nematic 
phases separated by first order phase transitions between 
them. However, to verify this assumption much larger 
system sizes would be required. 

Next, we considered the case 5 — 1, for which the 
ground state corresponds to a striped state with h = 1 
and the transition temperature is located aroundii T = 
0.4. In FiglHl we show the typical structure of the his- 
tograms for the energy and the order parameter observed 
for system sizes up to L = 48. While we do not see in 
this case a double peak structure in the energy per spin 
histogram (FiglSt) as in the 5 = 2 case, the broadening 
of the energy histogram around T = 0.4 and its skewed 
form at both sides of the transition, together with a slight 
double peak structure of the orientational parameter (al- 
though difficult to see, the numerical data in FiglHb show 
the presence of a minimum for T = 0.401 between 77 = 
and |?7| « 0.6) suggest a single first order phase transi- 
tion. Moreover, a direct inspection of the typical spin 
configuration associated with the different distributions 
indicate a direct phase transition from the striped to the 
tetragonal phase, without any trace of an intermediate 
nematic state (notice that the location of the maximum 
of P(|ry|) moves continuously towards one below the tran- 
sition point). This scenario will be confirmed by a direct 
thermodynamical analysis in the next section. However, 
based in the previously observed finite size effects for 
(5 = 2, the possible existence of a nematic phase in a nar- 
row range of temperatures for much larger system sizes 
cannot be excluded. 
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FIG. 7: Histograms for 5 = 2 and L = 64, for tempera- 
tures around the nematic-tetragonal transition, (a) Energy 
per spin; temperatures from left to right: T — 0.804, 0.805, 
0.807, 0.810 and 0.815. (b) Orientational order parameter; 
same temperatures as in (a), but from right to left. 



cific heat and the fourth order cumulant 

C{L,T) = j;^^ {{H') - (Hf) (11) 



V{L,T) = 1- ^ (12) 

To check the simulations results we also calculated the 
average energy per spin u(T) = — (H) /N and compared 
its derivative with respect to temperature with Ea. ljllll : 
both results coincide within the numerical error. We also 
analyzed the average of the absolute value of the orien- 
tational order parameter (|?7|) and the associated suscep- 
tibility 



III. PHASE TRANSITIONS 

We now turn our attention to the stripe-nematic and 
nematic-tetragonal phase transitions in the 6 = 2 case. 
To characterize them we analyzed the finite size scaHng 
behavior of the moments of the energy, namely, the spe- 



X,(L,T)^-((ry2)-(H)^) (13) 

All these quantities were obtained from the corre- 
sponding histograms for system sizes up to L = 64. In 
FigElwe show the behavior of the energy moments and in 
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FIG. 8: Histograms for 5 = 1 and L = 48. (a) Energy per 
spin; temperatures from left to right: T = 0.380, 0.390, 0.395, 
0.399, 0.401, 0.405and 0.410. (b) Orientational order param- 
eter; same temperatures as in (a), but from right to left. 
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FigElof the average order parameter and its associated 
susceptibility (the results for L = 64 are not shown for 
clarity) . 

We see that the specific heat shows two distinctive 
maxima at size-dependent pseudo-critical temperatures 
Tf (i) < T|(i), while the fourth order cumulant shows 
two distinctive minima at pseudo-critical temperatures 
Ti{L) < T2 (L), associated with the stripe-nematic and 
the nematic-tetragonal transitions respectively (notice 
that both the maximum at and the minimum at T2 
are almost unperceptive for L = 32). On the other hand, 
the orientational order parameter shows jumps in be- 
havior at each transition, but the associated suscepti- 
bility only shows a clear size dependent maximum at a 
size-dependent pseudo-critical temperature T2{L), cor- 
responding to the nematic-tetragonal transition. This is 
consistent with the fact the 90° rotational symmetry is 
broken at this phase transition. Although a second peak 
seems to emerge at lower temperatures for the largest size 
considered, it is very small and much more larger system 
sizes would be required to asses the presence of a size 
dependent peak associated with the stripe-nematic tran- 
sition. Prom FigEJwe see that the pseudo-critical tem- 
peratures of the specific heat and the fourth order cumu- 



lant scale as T^{L) - Ti -f Ai/L^, T^{L) - Ti + Bi/L^, 
T^{L) - + r|(L) - Ta + Bs/L^, with Ai > Bi 

and A2 > B2- This is the expected scaHng behavior 
at first order phase transitiona^fiiSl. The extrapolation 
to L ^ CO gives the values Ti = 0.772 ± 0.002 and 
T2 = 0.797 ± 0.005 showing that, although narrow, there 
is a well resolved range of temperatures in which the ne- 
matic phase exists in the thermodynamic limit. More- 
over, from Fig[T3we also see that T^{L) ^ T2 + C/L'^, 
where the extrapolation to L ^ 00 gives the same value 
for T2 as the previous calculation. This behavior is also 
consistent with a first order phase transition^i. 

However, when looking at Fig|2|a clear difference in fi- 
nite size scaling of the maxima of the specific heat and the 
minima of the fourth order cumulant between both tran- 
sitions appears: while the minimum of V at T2 rapidly 
saturates at a constant value when L — > oo, the minimum 
of V at Ti shows scaling behavior, approaching system- 
atically to 2/3. From Fig^jwe see that the maximum 
of the specific heat at increases as for system sizes 
up to L = 56, as expected in a first order transition, but 
it shows a little departure from this behavior for L = 64. 
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FIG. 10: Moments of the orientational order parameter 
Eq.(|2ll for 5 = 2 and different system sizes, (a) Average of 
the absolute value; (b) Associated susceptibility En. lllM|l . 
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FIG. 11: Finite size scaling of the pseudo-critical tempera- 
tures for the maxima of the specific heat (Tf and ) and the 
minima of the fourth order cumulant (Tf and T^) for 5 = 2 
and system sizes up to L = 64. The inset shows a zoom of 
the same curves around 
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However, such departure is probably due to systematic 
errors introduced by the presence of two phase transitions 
(see section InJ that cannot be resolved for small system 
sizes. The scaling of is similar to that of C(T|). 

On the other hand, the maximum of C at Tf saturates 
clearly in a finite value when the system size increases. 
It can also be shown that 2/3 - V{T^) ~ L'"^. This be- 
havior can be understood if we look at the energy per 
spin histogram at the temperature where both maxima 
have the same height. If we call estr{L) and CnemiL) the 
energies per spin at which those maxima occur, in a first 
order phase transition it is expected that they will scale 

as™ 6str{L) Ustr ^ L and &nem{^L^ Unera ^ L , 

Ustr and Unem being the specific internal energies of both 
phases (the striped and the nematic ones, in this case) 
in the thermodynamic limit. Then, the maximum of the 
specific heat is expected to scale in a two-dimensional 
system as^i C(Tf) ^ V- (e„em(L) - estr(L)) , while 
IIZ-V^TI) ~ (e„em(i)-e,t,(L)f + 0(L-2). In Figd 
we see that estr{L) and enem{L) show the expected 
scaling, but they converge to the same value, so that 
GnemiL) — CstriL) ~ and this explains the observed 
behaviors of the maximum of C and the minimum of 



FIG. 12: Finite size scaling of the pseudo-critical tempera- 
tures for the maxima of the susceptibility for 5 = 2 and system 
sizes up to 1/ = 64. 



V. This shows that the internal energy at the stripe- 
nematic transition becomes continuous in the thermo- 
dynamic limit. However, a thermodynamical singularity 
still develops in that limit. This can be seen by looking 
at the scaling of the free energy barrier between both 
phases, which can be calculated as^i 



AF ^ In 



Prr 



Prr 



(14) 



where Pmax = P{enem) = P{estr) is the value of the com- 
mon maximum of energy per spin histogram and Pmin is 
the value of the minimum between them. In Fig^l we 
see that AF ~ L in both phase transitions, which is 
the expected scaling in a first order phase transitional. 
In particular, in the case of the stripe-nematic transi- 
tion it means that, even when both states acquire the 
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FIG. 14: Finite size scaling of the energies per spin at which 
both maxima in the corresponding probability distribution 
occur, for the stripe-nematic phase transition for 5 = 2 and 
system sizes up to L = 64; for every value of L the tempera- 
ture is chosen such that P(e„em) — P(&str) 



same energy, there is a divergent free energy barrier be- 
tween them. Finally, we analyzed the difference between 
the values of the parameter 77 at which the maxima in 
the corresponding histogram occur in the stripe-nematic 
phase transition A7y(L) = r\str — Vnem- A finite size 
scaling analysis similar to the previous ones shows that 
Ari{L) ~ A77(oo) 4- D/L^, with Ar]{oo) = 0.075 ± 0.01. 
Hence, a finite jump in this parameter, even small, per- 
sists in the thermodynamic limit, consistently with the 
discontinuous jump in the stripe order parameter already 
observed in the previous section. 

A specific heat peak saturation has been observed nu- 
merically in different systems exhibiting a KT type phase 
transition, such as the 2D XY model^Si and generaliza- 
tions and the ID Ising model with 1/r^ ferromagnetic 
interaction922i21. This last model is particularly interest- 
ing, since it presents long range order in the low tem- 



perature phase (at variance with the XY models, which 
present only short range order) and the order parame- 
ter (magnetization) is discontinuous at the transitions^, 
in an analogous way to the model we are considering 
here. However, an important difference should be noted: 
in the above mentioned examples of the KT transition 
the extrapolation of the pseudo critical temperature at 
which the maximum of the specific heat is located ap- 
pears slightly above the critical temperature. 

Next we considered the tetragonal-stripe phase transi- 
tion in the case S — 1. As we have seen in section^ the 
equilibrium histograms suggest a weak first order phase 
transition. Since a finite size scaling approach would re- 
quire system sizes much larger than the allowed by the 
present computational capabilities, to confirm the sup- 
posed first order nature of the transition we considered 
another type of approach. 

Since we know the value of the ground state energy 
Ug and also the asymptotic infinite temperature inter- 
nal energy, we can obtain the free energy for different 
temperatures integrating the measured internal energies. 
We checked that the internal energy is independent of 
system size for 48 < L < 64, and we chose a iV = 48 x 48 
lattice in order to measure internal energy. The equilib- 
rium internal energy as a function of temperature for the 
stripe and tetragonal phases are shown in Fig. lfTfill . Using 
this energy curve we fit the stripe energy by a polyno- 
mial function UstriT) =Ug + aiT"^ + asT"-^ + a^T''^ and 
the tetragonal energy with a sum of hyperbolic functions 
Utetr{T) = &otanh(^) -I- 62tanh(^), whose parameters 
are shown in the caption. By means of the thermody- 
namical relation: 



= /3o/o(/3o) + / u{(})d(5, (15) 

where (3 = 1/T and Po = 1/To, we obtained the free 
energy per particle, /(/3) for the striped and tetragonal 
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FIG. 16: Internal energy per spin u versus T for 5 = 1 
and A'' = 48 X 48. Also shown best fits to the low tem- 
perature stripe and high temperature tetragonal regions, 
as described in the text. The best fit parameters for 
Ustr are ai = 1952.54; a2 = 10.8491; as = 10.8441; a4 = 
5.73016; as = -58.0638 and ae = 7.34061, and for Ut^tr 
are bo = -1.47926; bi = 0.0817682; 62 = 0.0559538 and 
63 ~ —3.05417. The ground state energy is Ug — —0.4677. 



phases (shown in Fig. Ill7|l l. The transition tempera- 
ture is obtained by imposing fstr{T) = ftetr{T), which 
gives Tm = 0.404. The entropy is obtained directly by 
the thermodynamical relation s = — ^ and the result is 
shown in the Fig. lll8ll . A weak first order transition is 
apparent, in agreement with the behavior of the energy 
and order parameter histograms of Fig|S| 

These results for the thermodynamic functions point 
to a unique phase transition for S = 1, from a high tem- 
perature tetragonal liquid phase to a low temperature 
striped phase. There is no trace of a third Ising nematic 
phase as for (5 = 2, at least for 48 < L < 64 system sizes. 

This conclusion is natural in this case where the equi- 
librium low temperature phase corresponds to stripes of 
width h = 1. Once defects begin to emerge due to ther- 
mal fiuctuations, disclinations appear naturally and the 
orientational order will rapidly decay , at variance with 
the expected behavior for wider stripes which will be 
more stable to the appearance of topological defects. 



FIG. 17: Free energy per spin / versus T for 5 = 1. The 
stripe and tetragonal free energies cross each other around 
Tm = 0.404. The functions obtained from equation Ijl5ll 

are = + ^r"= + ^T"^ + ^T''^ and ft.u- = 

l^rin(cosh^) + ^rin(cosh^) - rin(2), where /3o = for 
the tetragonal and Po = oo for the striped phase. 
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FIG. 18: Entropy per spin s versus T for 5 = 1. The slight 
discontinuity at the transition temperature is consequence of 
a weak first order phase transition. 



IV. METASTABILITY 

The presence of diverging free energy barriers between 
the different phases found in the previous section for 
S = 2 lead us to consider the possible existence of dy- 
namical metastable (i.e., quasi-stationary) states. To 
this end, we performed cooling and heating numerical 
calculations according to the following protocol. We first 
let the system equilibrate at a temperature high enough 



to ensure that it is in the tetragonal Hquid phase. We 
then cooled the system at a fixed cooling rate r, that is, 
the temperature was decreased as T{t) = Tq — r t, where 
To is the initial temperature and the time t is measured 
in MCS. The temperature was reduced in every MC run 
down to a value well below the range associated with the 
different phase transitions, while recording at different 
points the energy and the orientational order parameter 
\t]\ of the system. Every curve was then averaged over 
different initial conditions and different sequences of the 



random noise. For every system size we calculated those 
curves for decreasing cooling rates until they became in- 
dependent of the cooHng rate. In this way we simulated 
a process of quasi-static cooling. Once we obtained the 
quasi-static cooling curves we performed a quasi-static 
heating from zero temperature up to high temperatures, 
starting from the ground state and using the same pro- 
tocol and the same rate r. 

In fig^lwe see an example of the quasi-static cooling- 
heating curves obtained for L = 48 with a cooHng rate 
r = 10~^; the results are compared with the equihbrium 
curves obtained in the previous section. Similar results 
were obtained for system sizes up to L = 64. We do 
not observe super cooled tetragonal Hquid states below 
T2, but this is probably due to the relatively small val- 
ues of the associated free energy barriers for the system 
sizes here considered. On the other hand, we observe a 
strong metastability associated with the stripe-nematic 
phase transition, which is consistent with the observed 
large free energy barriers. In particular, we see that the 
supercooled nematic state is observed down to temper- 
atures well below Ti. To verify the quasi-equihbrium 
nature of the metastable nematic states, we calculated 
the two-times autocorrelation function 



(16) 



after a quasi-static cooHng to different temperatures 
T < Ti, down to T = 0.5, for different pairs of values 
ti <t2. In all the cases we found that C{ti,t2) depends 
only on the difference t2 — ti (as expected for a quasi- 
stationary process), for time scales up to t = 10^ MCS. 
For larger time scales the nematic phase finally decay into 
the striped one, apparently through a nucleation process. 
These results are not shown here, but a complete descrip- 
tion of those calculations will be published shortlji2£. 



V. DISCUSSION 

The main contribution of this paper are the several 
evidences of the existence of nematic phases in an ul- 
trathin magnetic film model, at least for some range of 
values of the exchange to dipolar intensities ratio. These 
results are in agreement with one of the two possible 
scenarios of the critical behavior of those systems pre- 
dicted by the continuum approximation of Abanov and 
co-workerafi. Let us discuss first the stripe-nematic phase 
transition. Although the finite size scaling behavior of 
different quantities around that transition agree with 
that expected in a first order transition, both the fact 
that the energy becomes continuous and the saturation 
of the specific heat maximum in the thermodynamic limit 
are unusual in a simple first order transition, suggesting a 
more complex phenomenology. In fact, those properties 
strongly resemble those observed in the one dimen- 
sional Ising model transition point, namely, continuous 
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FIG. 19: Cooling-heating curves for 5 = 2, L = 48 and 
r = 10~®. The insets show the same curves for a wider range 
of temperatures, (a) Average energy per spin (b) Average of 
the absolute value of the orientational order parameter. 



energy, discontinuous order parameter and saturation of 
the specific heat maximum. These analogies suggest that 
some KT type mechanism, probably associated with the 
unbinding of dislocations, could be responsible of the ob- 
served phenomenology, in agreement with Abanov and 
coworkers prediction. The observed first-order-like finite 
size scaling behavior appears to be related with the dis- 
continuous change of the orientational order parameter at 
the transition point in the thermodynamic limit, which 
suggests the existence of a finite density of dislocation 
pairs. Indeed, we observed the presence of an increas- 
ing number of dislocations pairs (bridges) in the striped 
phase at r = Ti but, up to L = 72 system sizes we didn't 
find any evidence of QLRO (algebraic decaying correla- 
tions). One important consequence of such first-order- 
like phenomenology is the existence of diverging free en- 
ergy barriers at that point, which has an important phys- 
ical implication, namely, the possibility to obtain highly 
stable supercooled nematic states through cooling. This 
opens the possibility of having complex slow dynamical 
behaviors after a sudden quench, such as those observed 
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in supercooled molecular liquids. We are investigating 
this problems and the results will be presented in a forth- 
coming publicationp2&. 

Concerning the nematic-tetragonal phase transition, 
our results suggest a much more complex behavior than 
the expected from Abanov et al. results, who conjec- 
tured a second order phase transition. Although our re- 
sults are inconclusive due to the presence of strong finite 
size effects, they suggest a complex phase diagram with 
multiple nematic phases, separated by first order tran- 
sition fines, similar to that encountered by Grousson et 
al2i in a related model, the 3D Coulomb frustrated Ising 
ferromagnet. However, the departure from the expected 
first order finite size scafing observed for L = 64 in the 
maximum of the specific heat (see Fig lT^jl and the sus- 
ceptibility, could be an indication of a behavior similar 
to that observed in the stripe-nematic transition. Hence, 
the possibility of a nematic-tetragonal KT phase tran- 
sition, associated with a disclination unbinding mecha- 
nism, cannot be excluded. A similar scenario appears in a 
related problem, namely, the two-dimensional melting^S,. 
On the other hand, the Hartree approximation of the 
Landau-Ginzburg version of the present model predicts 
a fiuctuation induced first order phase transitions^ for 
any value of 5. 

The case of small exchange/dipolar ratio is simpler in 
the sense that only one transition is observed from a 
stripe phase directly to a tetragonal liquid phase. We 
obtained evidence that this transition is a weak first or- 
der one. 

Although we have analyzed the critical behavior of this 



system for relatively small values of the exchange/dipolar 
ratio {6=1 and i5 = 2), the present results may help to 
understand the critical behavior for larger values. Booth 
and coworkers^ analyzed the behavior of the specific heat 
and the orientational order parameter for (5 = 3 and 
S = 4.45 and system sizes up to L = 64. While their re- 
sults suggest a continuous stripe-tetragonal phase transi- 
tion, numerical results based on time series analysis from 
Casartelli and coworkersiS, for the same parameters and 
system sizes, indicate that the transition is first order. 
For those values of 6 the stable low temperature phase 
is composed by stripes with a larger width {h = 4), so 
much more larger system sizes would be required to allow 
the appearance of a large density of topological defects 
(dislocations and disclinations of the stripes). Therefore, 
for larger system sizes one could expect a scenario simi- 
lar to that observed for 6 = 2, namely, the appearance of 
nematic order in a narrow range of temperatures. 

The experimental verification of the different transi- 
tions and phases suggested by theory and simulations is 
an important challenge which can be reached in the near 
future due to the emergence of new observation tech- 
niques. 

Fruitful discussions with A. Cavagna and T. Grigera 
are acknowledged. We thank very useful suggestions and 
comments from anonymous referees. This work was par- 
tially supported by grants from CONICET (Argentina), 
Agenda Cordoba Ciencia (Argentina), SeCyT, Univer- 
sidad Nacional de Cordoba (Argentina), CNPq (Brazil) 
and ICTP grant NET-61 (Italy). 



* Electronic address: cannasQfamaf.imc.er lii .arl 
t Electronic address: michelon(§lifi.uni camp.br| 

* Electronic address: stariolo@if.ufrgs. br| 

^ Electronic address: tamaritQfamaf.imc.e rin.arl 
^ Member of CONICET, Argentina 

** Present address: Institute de Fisica Gleb Whataghin, UNI- 

CAMP, Campinas, Brazil 
^ R. Allenspach, M. Stampanoni, and A. Bischof, Phys. Rev. 

Lett. 65, 3344 (1990). 
^ A. Vaterlaus, C. Stamm, U. Maier, M. G. Pini, P. Politi, 

and D. Pescia, Phys. Rev. Lett. 84, 2247 (2000). 
^ K. De'Bell, A. B. Maclsaac, and J. P. Whitehead, Rev. 

Mod. Phys. 72, 225 (2000). 

P. Politi, Comments Cond. Matter Phys. 18, 191 (1998). 
^ I. Booth, A. B. Maclsaac, J. P. Whitehead, and K. De'Bell, 

Phys. Rev. Lett. 75, 950 (1995). 
" A. Abanov, V. Kalatsky, V. L. Pokrovsky, and W. M. 

Saslow, Phys. Rev. B 51, 1023 (1995). 

L. C. Sampaio M. P. deAlbuquerque and F. S. deMenezes, 

Phys. Rev. B 54, 6465 (1996). 

A. D. Stoycheva and S. J. Singer, Phys. Rev. Lett. 84, 4657 
(2000). 

^ J. H. Toloza, F. A. Tamarit, and S. A. Cannas, Phys. Rev. 

B 58, R8885 (1998). 
" D. A. Stariolo and S. A. Cannas, Phys. Rev. B 60, 3013 
(1999). 



P. M. Gleiser, F. A. Tamarit, and S. A. Cannas, Physica 
D 168-169, 73 (2002). 

P. M. Gleiser, F. A. Tamarit, S. A. Cannas, and M. A. 
Montemurro, Phys. Rev. B 68, 134401 (2003). 
S. A. Cannas, D. A. Stariolo, and F. A. Tamarit, Phys. 
Rev. B 69, 092409 (2004). 

M. Casartelli, L. Dall'Asta, E. Rastelli, and S. Regina, J. 
Phys. A: Math. Gen. 37, 11731 (2004). 
O. Portmann, A. Vaterlaus, and D. Pescia, Nature 422, 
701 (2003). 

" Y. Z. Wu, C. Won, A. SchoU, A. Doran, H. W. Zhao, X. F. 
Jin, and Z. Q. Qiu, Phys. Rev. Lett. 93, 117205 (2004). 
O. Portmann, A. Vaterlaus, and D. Pescia, Phys. Rev. 
Lett. 96, 047212 (2006). 

It is worth noting that the definition of the Hamiltonian 
in Referencea^Aifii is slightly different from ours. While in 
those papers the dipolar term contains a sum over all pairs 
of spins, here we consider the sum over every pair of spins 
just once. This leads to the equivalence 5 = J/2, J being 
the exchange parameter in the above references. Since the 
dipolar parameter also fixes the energy units in our work, 
there is also a factor 1/2 between the critical temperatures 
obtained in those works and ours. 

" A. B. Maclsaac, J. P. Whitehead, M. C. Robinson, and 
K. De'Bell, Phys. Rev. B 51, 16033 (1995). 

^° M. S. S. Challa, D. P. Landau, and K. Binder, Phys. Rev. 



13 



B 34, 1841 (1986). 

J. Lee and J. M. Kosterlitz, Phys. Rev. B 43, 3265 (1991). 
R. Gupta and C. F. Baillie, Phys. Rev. B 45, 2883 (1992). 
E. Luijten and H. Messingfeld, Phys. Rev. Lett. 86, 5305 
(2001). 

J. Z. Imbrie and C. M. Newman, Commun. Math. Phys. 
118, 303 (1988). 

M. Ainzenman, J. T. Chayes, L. Chayes, and C. M. New- 



man, ,J. Stat. Phys. 50, 1 (1988). 

S. A. Cannas, D. A. Stariolo and F. A. Tamarit, to be 
published. 

M. Grousson, G. Tarjus, and P. Viot, Phys. Rev. E 64, 
036109 (2001). 

K. J. Strandburg, Rev. Mod. Phys 60, 161 (1988). 



